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Laser induced ultrafast demagnetization is the process whereby the magnetic moment of a ferro¬ 
magnetic material is seen to drop signihcantly on a timescale of 10 — 100s of femtoseconds due to 
the application of a strong laser pulse. If this phenomenon can be harnessed for future technology, 
it offers the possibility for devices operating at speeds several orders of magnitude faster than at 
present. A key component to successful transfer of such a process to technology is the controllability 
of the process, i.e. that it can be tuned in order to overcome the practical and physical limitations 
imposed on the system. In this paper, we demonstrate that the spin-orbit mediated form of ultrafast 
demagnetization recently investigated arXiv: 1406.6607 by ab-initio time-dependent density func¬ 
tional theory (TDDFT) can be controlled. To do so we use quantum optimal control theory (OCT) 
to couple our TDDDT simulations to the optimization machinery of OCT. We show that a laser 
pulse can be found which maximizes the loss of moment within a given time interval while subject to 
several practical and physical constraints. Furthermore we also include a constraint on the fluence of 
the laser pulses and hnd the optimal pulse that combines significant demagnetization with a desire 
for less powerful pulses. These calculations demonstrate optimal control is possible for spin-orbit 


mediated ultrafast demagnetization and lay the 
which can incorporate even more constraints. 

INTRODUCTION 

Control of quantum dynamics using tailored laser 
pulses is a long standing goal of modern physics [T]-[7] 
as it opens up a whole new world of possibilities for 
future technologies. Faster, smaller, and more efficient 
devices could be constructed if we could master con¬ 
trol over the charge and spin dynamics of electrons on 
the nanoscale [8]. However precisely at these very short 
length and time scales, quantum effects are strong, which 
makes it difficult to exert this control. With the advent [9] 
of laser pulse shapers that can tailor the laser field to a 
given shape, there was now a tool that could be used for 
control of quantum dynamics. The challenge is finding 
the shape of the laser pulse that produces the desired 
dynamics. 

Optimal control theory (OCT) is a method 
developed [TOl [H] in both Mathematics and Engi¬ 
neering to solve the problem of finding a particular 
control variable that gives a desired outcome. In our 
case, we will search for the electric field e(t) of a laser 
pulse to control the properties of our system. In general 
OCT works by creating a target functional of the control 
field calculated from simulation of the system. Then any 
constraints on the system are incorporated using penalty 
functionals, before extremizing the total functional to 
find the optimal field. OCT can be extended to the 
realm of quantum mechanics by constructing the target 
functional using observables given by the time-dependent 
schrodinger equation (TDSE). 

Eor more than a handful of electrons, propagating the 
TDSE is a computationally intractable problem due to 


foundation for future optimizations/simulations 


the coulomb interaction between electrons and an alter¬ 
native approach must be used. Time-dependent density 
functional theory (TDDET) is one such approach, which 
works by mapping the problem to a non-interacting 
system[T2], referred to as the Kohn-Sham (KS) system. 
This system is defined such that propagating electrons 
in this system will reproduce the same time dependent 
density (the probability to find an electron at any given 
point) as propagating in the exact system using the 
TDSE. As the KS system is non-interacting, the prob¬ 
lem is now computationally tractable. Although, in prin¬ 
ciple, this mapping is exact, in practice an approxima¬ 
tion must be used. The most commonly used approxi¬ 
mation, adiabatic local density approximation (ALDA), 
has been shown to successfully predict absorption spec¬ 
tra of a large range of atoms, molecules, and solids [13]- 
[I5]. Thus TDDET is an outstanding candidate to couple 
to OCT [16] in order to predict laser pulses for control 
of quantum dynamics, and has been used successfully 
for control of charge transfer[T7]. HHG[20|, strong-field 
ionization[T8l [19], bond-breaking [21], among others. 

Laser-induced ultrafast demagnetization was first ob¬ 
served in the mid 1990s, whereby a strong femtosecond 
laser pulse caused a significant loss of the magnetic mo¬ 
ment of a thin film of Ni in a time less than Ids [22]. 
Since then, this phenomena has been the subject of much 
experimental [23H33] and theoretical [34H37] endeavor and 
several mechanisms have been proposed to explain the 
demagnetization. In Ref. [371 ab-initio TDDET simula¬ 
tions were performed to investigate the demagnetization 
and found that when spin-orbit interaction was included 
in the system Hamiltonian, a loss of moment was ob- 
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served for very short (5fs), very intense (1 x 10^^ W/cni2) 
laser pulses. It is this system we wish to control by vary¬ 
ing the intensity and frequency of the laser pulse, subject 
to several practical constraints, in order to maximize the 
total loss of moment. To do so we utilize the framework 
developed in Refs. |T6l [2TJ and [38l which combines OCT 
with quantum simulations of spin dynamics. 

BACKGROUND AND METHODS 


The KS potential, i;s(r, t), consists of three pieces: 

Vs{r,t) = (r, t) + i;h (r, t) + i;xc (r, t) (6) 

where is the usual Hartree potential of the in¬ 

stantaneous density, and r’xc(i*,t) is the exchange corre¬ 
lation (XC) potential and is a functional of the density 
at all previous times, the interacting initial state, and 
the non-interacting initial KS state. In practice, it must 
be approximated, with the most common approximation 
being the ALDA: 


We begin by briefly reviewing TDDFT and OCT, a 
more thorough discussion can be found here [16]. Then 
we review the results of Ref. [37] that showed ultrafast de¬ 
magnetization in bulk ferromagnets (Fe,Co,Ni) for short, 
strong, laser pulses. 


TDDFT 

The electronic density is defined as 

n{r,t) = N J dr 2 . ..drN 4'*(r,r2, 

x^'(r,r 2 ,...,rjv,t) (1) 

where N is the total number of electrons, r is the spacial 
coordinate, t is the time, and T is the wavefunction of 
the TDSE: 


t>xcM(r,t) = 


de 


unif 

XC 


dn 


n=n(r,t) 


( 7 ) 


which uses just the instantaneous density inputed into 
the ground-state DFT LDA XC functional and {n) 
is the XC energy density of the uniform electron gas. The 
initial KS state is typically the ground-state found from 
a DFT calculation. 

From this starting point, TDDFT has been extended to 
include non-collinear magnetism and magnetic fields [40]. 
For this case, we have a non-interacting Pauli KS 
Hamiltonian!^ which is used to propagate 2 compo¬ 
nent spinors, from which the density and magnetization 
density exactly replicate those of the interacting system. 
This is the formulation we will use for our simulations. 
The magnetization density operator may be written as: 


for Hamiltonian: 


II 

(2) 

m(r) = S h(r) 

(8) 

H — T ^ Wxt + Vee 

(3) 

where h(r) is the density operator 
component spinors propagated in our 

and in the two- 
calculations, S = 


composed of the kinetic energy, T, the electron-electron 
interaction, and the external potential, Kxt? which 
includes both the electron-nuclear interaction and the 
electric fields of any laser pulses. We use atomic units 
throughout unless otherwise stated. TDDFT is founded 
upon the Runge-Gross theorem [T2] which proves a 1 — 1 
correspondence between the time-dependent density and 
the time-dependent external potential (up to a time- 
dependent constant) for any electron-electron interac¬ 
tion. Hence all observables of the system are, in prin¬ 
ciple, unique functionals of the density. In particular, 
a non-interacting KS system [39] can be defined with a 
unique KS potential that reproduces the time-dependent 
density of the interacting system and thus predicts all ob¬ 
servables of the true system without requiring the costly 
propagation of Eq. The TDKS equation is: 

~+Vs{r,t) 


d 


( 4 ) 


with the total density given by 

N 

i=i 


( 5 ) 


gj^cr where <j^, cr;^} are the familiar Pauli spin ma¬ 
trices and g is the electronic gyromagnetic ratio. Eor 
periodic boundary conditions, the total moment is then 


M(t) = f d^rm(r, t) 

Jn 


(9) 


where D is a single unit cell. The KS Hamiltonian for 
our simulations is: 

Hs{t) = ^(p+ ^Aext(^))^ + t’s(r,t) 

+ is • Bs(r, t) + ^S- {Vvs{r,t) x p) (10) 

where p is the momentum operator, S is the vector spin 
operator, and c is the speed of light. The laser pulse 
electric field is written as a vector potential, Aext (t) in 
the velocity gauge as it allows Bloch’s theorem to be uti¬ 
lized. The KS magnetic field is written as Bs(r,t) = 
Bext(^) + BxcCf*,^ where Bext(^) is the magnetic field 
of the applied electromagnetic field and Bxc(f*,^) is the 
XC magnetic field. The ALDA can be extended to Bxc 
using the LDA rotation method of Kubler[42]. The final 
term of Eq. (10) is the spin-orbit coupling (SOC) term. 
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which can be thought of as the interaction between the 
spin of an electron and the effective magnetic field caused 
by relativistic motion thought a scalar potential. In a 
centrosymmetric potential, this term reduced to the well 
known L • S coupling. Propagation with Hamiltonian Eq. 
(10) is implemented in ELK [43], an all-electron electronic 
structure code, which was also used for all ground state 
and time-dependent calculations. 


Optimal Control Theory 

The central equation of OCT is the target functional 
G[u]: 


G[u]=Gmu],u] = Ji[u] + J 2 [u] (11) 

where u is the control field and 4^ [u] contains the informa¬ 
tion on how the system responds to the control field. In 
Quantum OCT (QOCT), is then the wavefunction, 
which is a functional of the control field via the TDSE 
and from which any system observables to be controlled 
may be calculated. The target functional is generally 
separated into two pieces, Ji[u] which contains informa¬ 
tion on the desired dynamics and J 2 [u] which is a penalty 
function in order to satisfy any constraints on the system 
or control field. The magnitude of the penalty functional 
is determined by how strongly a constraint must be sat¬ 
isfied. 

Once a relevant target functional has been constructed, 
the goal of OCT is to extremize it and thus find the 
optimal control field u to best satisfy the balance be¬ 
tween desired dynamics and the constraints. There are 
many choices for the algorithm to perform this opti¬ 
mization, some are general, such as the Nelder-Mead [44] 
or NEWOAU[45] algorithms, while some are developed 
for specific types of problem, e.g. in QOCT the ZBR 
scheme [46] adds a time dependent auxiliary wavefunc¬ 
tion, which is also propagated in time and the overlap 
with the true wavefunction used to construct the control 
field. 

Eor our system, we wish to maximize the loss of mag¬ 
netic moment in a given time interval [0,T] while includ¬ 
ing practical and physical constraints on the type of laser 
pulse. Thus e(t), the electric field of the laser pulse is the 
control field and 

Me] = (^[e](T)|M,|^[e](T)) = M,(T) (12) 

is the target functional to be minimized, i.e. if we choose 
the initial magnetization Mz{0) of the ferromagnet to 
be along the 2 )-axis, then minimizing Mz{T)/Mz{0) will 
maximize the loss of moment. 

The constraints on the electric field are that the pulses 
satisfy Maxwell’s equations (details below) and only cer¬ 
tain frequencies are used to construct the pulse. The sec¬ 
ond constraint is of practical nature, as experimentally. 
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FIG. 1. Upper panel: The profile of the laser vector poten¬ 
tial (labelled A field). Lower Panel: The total moment per 
atom as a function of time for several different peak intensities 
(given in W/cm2). [Figure reproduced from Ref. [37l ] 


pulses containing arbitrary frequencies cannot be con¬ 
structed and often access to a single frequency (or mul¬ 
tiples thereof) is only available. From Maxwell’s equa¬ 
tions, the following constraints on the electric field must 
be physically satisfied: 


/ 


dte{t) = 0 


6(0) = 0 = e(T) 


(13) 

(14) 


Following Ref. EH we can satisfy all these constraints by 
writing the electric field 


e{t) = y] enCOs(aJnt) (15) 

n=l 


where is the number of frequencies to be used and 

in are the coefficients to be optimized. It can be seen 
that this choice automatically satisfies the constraints 
from Maxwell’s equations. The frequencies used for our 
demonstration are 


UJn 


27vn 

~T~ 


(16) 


Ultrafast Demagnetization 

In Fig. the results of TDDFT simulations [37] for 
several laser pulses with differing peak intensities but the 
same profile are shown. A loss of the total magnetic 
moment was observed in all cases, with the fraction of 
moment lost dependent on the field intensity. It was also 
shown in Ref. EH that the fraction loss depends on the 
frequency of the applied laser pulses. Hence the system 
was a strong candidate for optimal control. The purpose 
of this paper is test that hypothesis. 

It should be pointed out that for less-intense longer- 
duration pulses and for more realistic system geometries, 
the required peak intensity and fluence can be reduced by 
several orders of magnitude [47], however the underlying 
mechanism of demagnetization is the same. Thus we can 
demonstrate control of this process by focusing on the 
short strong laser pulses. 
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all quantities in atomic units 



FIG. 2. Upper Panels: The vector potentials for initialization 
laser pulses. Lower Panels: The dynamics of the total moment 
in the ; 2 ;-direction for each pulse. 


RESULTS 


To demonstrate optimal control of the ultrafast demag¬ 
netization in Ni, we will attempt to maximize the loss of 
moment after T = 600 au = 14.4 fs. We will optimize a 


pulse of the form given by Eq. (15) using = 4 different 


laser frequencies defined by Eq. (16). All calculations are 


performed with a time step of 0.1 au and 8x8x8 /c-points. 
Eor the optimization we choose to use the gradient-free 
Nelder-Mead simplex algorithm. To initialize the Nelder- 
Mead algorithm, +1 starting points are required, it is 
instructive to examine these before moving to the results 
of the optimization. 


Random Initial Pulses 


To initialize our calculation, we construct 5 different 


pulses where the coefficients of Eq. (15) are chosen at 


random in a suitable range. This range is chosen such 
that the peak intensity is similar to the demagnetization 
pulses observed in Ref. |371 The pulses may be seen in 
the upper panel of Eig. note that this is the vector 
potential which can be calculated from the electric field 


via 


A(t) = —c j dt' e{t') (17) 

The dynamics of (t) are shown in the lower panel and 
it can be seen that all pulses display demagnetization. 
If we look at the final time, the average loss of moment 
is approximately 20%. If the optimal control is success¬ 
ful, then this percentage loss should be significantly in¬ 
creased. 


Maximize Demagnetization 

Erom the initial pulses, the Nelder-Mead algorithm 
then calculates a new set of coefficients from a simple 
set of rules and then tests how this affects the target 
functional by performing a TDDET simulation with the 



FIG. 3. The fraction of magnetic moment loss for each itera¬ 
tion of the Nelder-Mead optimization. 


laser field given by these coefficients. It then iterates this 
procedure and traverses the multidimensional parameter 
space, searching for the optimal set of coefficients. 

In Fig. we plot the ratio of the final moment after 
time T to the initial moment for each of the iterations. 
Although individual iterations can worsen the loss per¬ 
centage, there is a clear downward trend as better and 
better pulses are found during the search, indicating that 
the optimal control is working. Each set of coefficients 
is a point in the parameter space, at each iteration, the 
Nelder-Mead algorithm reflects the worst point through 
the center of mass of the other points. Depending on 
whether this new point improves upon the next worst 
point, the algorithm can expand or contract in this direc¬ 
tion, otherwise it can reduce all points towards the best 
point. This explains why individual points may worsen 
the ratio M^(T)/M; 2 ( 0 ). 

If we look at the result after approximately 30 itera¬ 
tions, the best pulse the optimal control procedure has 
found causes a 40% loss of moment. This is twice as good 
as the random initial pulses used to start the algorithm 
and is a clear demonstration that the moment can be 
successfully controlled using OCT. In Eig. [^we show the 
electric field of this best pulse and also the magnetization 
dynamics, compared to the initial pulses. Examining the 
pulse shape compared to the initial pulses, there is no ob¬ 
vious reason why one leads to a larger demagnetization. 
This is the power of QOCT to find such pulses. We also 
observe that the magnetization dynamics is a highly non¬ 
linear process, in particular for these short pulses, again 
demonstrating the need for QOCT in order to control the 
moment. 
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FIG. 4. Upper Panel: The optimal laser field found after 30 
iteration of the optimal control algorithm. Lower Panel: The 
magnetization dynamics of for this pulse compared to the 
5 initialization pulses. 


FIG. 5. The OGT target functional, J[u] = 
Mz[u]{T)/MzifS) + 0.05 X F[u] for the electric held u = e{t) 
at each iteration of the algorithm, where F is the huence. 


Fluence Constraint 


For many practical reasons, the fluence of the applied 
pulse should be constrained. The simplest reason is sim¬ 
ply efficiency, i.e. using a a pulse with lower energy 
to achieve the same dynamics as a higher energy pulse. 
Other reasons include surface damage to the material due 
to high fluence pulses, heating of the sample (and prob¬ 
lems associated with cooling it), or physical restrictions 
on the laser itself preventing production of high fluence 
pulses. All of these present significant problems to future 
technological application, hence we wish to demonstrate 
how a fluence constraint can be incorporated into our 
calculations. 

If we add to Eq. 0 , the constraint 

J2[e] = a f dt e^{t) (18) 

Jo 


which is proportional to the laser fluence. The free pa¬ 
rameter a determines how strong the constraint is, for 
this calculation we choose a = 0.05. This parameter was 
based on examining the results of the previous optimiza¬ 
tion and choosing a to favor a lower fluence while still 
maintaining significant demagnetization in the set. 

In Fig. we show the value of the total target func¬ 
tional, Eq. 0, for each iteration of the optimization 
algorithm. Unlike the previous case, we cannot attach 
a physical meaning to the target functional, so the ac¬ 
tual value is not significant, only the trend. Eurthermore, 
when choosing a, it was clear that the parameter space is 
a more complicated environment than the previous case, 
as a pulse could have the same value of Eq. (11) by 
either increasing the demagnetization or decreasing the 


fluence. Due to computational constraints, we stopped 
the optimization after 12 iteration, although this was suf¬ 
ficient to see the trend and demonstrate optimal control. 
We initialize the algorithm using the 5 best pulses found 
by calculating the target functional for all the previous 
pulses. These are not shown in Eig. but were part of 
the optimization search. By using these to start the opti¬ 
mization, we save a large amount of computational time 
as opposed to using random pulses (although in general 
this may not be the case if the constraints significantly 
change the parameter landscape). 

To see the power of this optimization, in Eig. |^we plot 
the electric fields and the dynamics of (t) for two dif¬ 
ferent pulses. These correspond to the best point of Eig. 
I^and a reference pulse corresponding to the G = 1.1021 
point. While this point is not the worst point, it was 
chosen as the final moment is similar to the best point, 
—0.537 and —0.555 respectively. Thus, the magnetic mo¬ 
ments at time T are very close, however the fluences are 
very different. If we insert the pulses shown in Eig. 
into the fluence formula given in Eq. we find values 
of 6.437 a.u for the reference pulse and 2.437 a.u. for the 
best pulse. Therefore by using OCT we have reduced the 
required fluence by ~ 60%. 


CONCLUSIONS 

In summary we have successfully demonstrated that 
optimal control of spin-orbit mediated ultrafast demag¬ 
netization is possible. Eor a short time interval, we 
showed that the loss of moment can be at least dou¬ 
bled (compared to randomly chosen typical pulses) for 
a system where the available laser frequencies (used to 
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FIG. 6. The electric field E(t) (upper panel) and total mag¬ 
netic moment (lower panel) for the reference pulse (solid 
black line) and the best pulse (dashed red) found during the 
fluence constraint optimization. 

tailor the laser pulse) are constrained. Furthermore we 
extended the control problem to include a constraint on 
the laser fluence and demonstrated that QOCT could 
successfully find a pulse that balances the fluence and 
demagnetization requirements. Compared to a reference 
pulse, this optimal pulse produces almost identical mag¬ 
netization dynamics, while reducing the fluence by over 
a factor of 2. Control of the system is of upmost impor¬ 
tance for future technological application (for example, in 
spintronics), where the desired dynamics and constraints 
are dictated by practical concerns. Any physical phe¬ 
nomenon must be robust to these concerns, and as we 
have demonstrated, this form of ultrafast demagnetiza¬ 
tion meets this criteria. Simulation and QOCT of more 
complicated scenarios, such as longer pulse durations or 
further constraints on the fluence, intensity, and robust¬ 
ness of the demagnetization, can be build upon this foun¬ 
dation. 
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